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We consider the non-equilibrium physics induced by joining together two tight binding fermionic 
chains to form a single chain. Before being joined, each chain is in a many-fermion ground state. 

The fillings (densities) in the two chains might be the same or different. We present a number 
of exact results for the correlation functions in the non-interacting case. We present a short-time 
expansion, which can sometimes be fully resummed, and which reproduces the so-called ‘light cone’ 
effect or wavefront behavior of the correlators. For large times, we show how all interesting physical 
regimes may be obtained by stationary phase approximation techniques. In particular, we derive 
semiclassical formulas in the case when both time and positions are large, and show that these are 
exact in the thermodynamic limit. We present subleading corrections to the large-time behavior, 
including the corrections near the edges of the wavefront. We also provide results for the return 
probability or Loschmidt echo. In the maximally inhomogeneous limit, we prove that it is exactly 
gaussian at all times. The effects of interactions on the Loschmidt echo are also discussed. 
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Introduction - Local quantum quenches are a par¬ 
ticularly neat setup to address theoretical questions 
about non-equilibrium current-carrying stationary states 
as well as transport properties of many-body isolated 
quantum systems. The characterization of the long-time 
behavior of local correlation functions unveils universal 
features of the quantum dynamics m and paves the 
way to the construction of effective field theories capable 
of capturing them. Analytic results are relevant as bench¬ 
marks for cold atom experiments that very recently EHH] 
started to investigate particle and energy transport under 
a unitary dynamics. 

In this work, we study a quench which, although in¬ 
jecting an extensive amount of energy into the system, 
has mostly local effects. We take two uniformly filled long 
tight binding chains in their respective ground-states, but 
with a different number of fermions. At time t = 0 the 
two edges are connected so that it turns into a single 
chain 

^ L/2-1 

^ = “2 X ■ (1) 

j=-L/2+l 

The hopping and interactions between sites j = 0 and 
j = 1 are initially absent; the initial state |/o) = IV’;) IV’r) 
is the tensor product of the ground states (with specihed 
occupancies) of the two decoupled chains. For unequal 
fillings, one expects some particle current at time t > 0. 
We denote by kp (kp) the Fermi momenta on the left 

(right), so that the particle number is 
left (right). For simplicity we focus on the symmetric case 
kp + kp = TT, but extension to other values is straight¬ 
forward. It is useful to keep two limits in mind. When 
the fillings are equal kJ'p = k'p = tt/2 there is no particle 
current. This particular quench was studied in [THSIE], 
using low-energy held theory, and belongs to the class 


of Fermi-edge problems [iniiii]. The other simple limit 
is k^p = TT, with a domain wall (DW) initial state [T^ 
El iV’o) = n^<o 4 |0), where |0) is the fermion vacuum. 
For intermediate Filing a non equilibrium steady state 
(NESS) develops in the middle of the chain, in which 
correlations are that of a Fermi sea k G [—kp, kp] with 
shifted momenta. This central region may also be under¬ 
stood using bosonization m- 

The motivation for the present study is two-fold. First, 
our initial state is more rehned than previously studied 
“projected Fermi seas” [HIIIIS], and fully accounts for 
the boundary effects near the junction. Second, we derive 
a number of exact results after the quench, in many inter¬ 
esting physical regimes. These include the long-time limit 
of correlations at hnite x, the boosted Fermi sea regime, 
where we are also able to derive the leading corrections. 
One of our main results concerns the regime at large po¬ 
sitions x,y and time t, with hnite ratios xlt,yjt. It is 
known that such regimes are qualitatively well described 
by semiclassical arguments. One complication in our case 
is that we are dealing with a discrete lattice model with 
non polynomial dispersion, where usual techniques based 
on the use of Wigner functions (see e.g [TTHI^ ) are not 
easy to generalize. We show nevertheless, by a careful 
use of stationary phase arguments in presence of singu¬ 
larities, that such a picture becomes exact in the thermo¬ 
dynamic limit, for |x —y| <C t. The present work provides 
a simple unihed picture of all physical regimes, in terms 
of initial state correlations and single particle energies. 
We also consider the Loschmidt echo, and combine ana¬ 
lytical and numerical techniques to treat the short- and 
large-time behaviors. The Loschmidt echo is also consid¬ 
ered in the presence of nearest-neighbor interactions. 

Power series expansion — We start by deriving an ex¬ 
act power series representation for the two point func- 
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tion. All higher order correlations may be obtained from 
Wick’s theorem. Consider two fermionic quadratic forms 
^ and H = where A 

and H are Lx L matrices. As is well known, see e.g. pO] . 
quadratic forms in Fock spaces are a representation of the 
Lie algebra gl{L), namely [A,H] = 

Interpreting H as the final Hamiltonian and A as an 
observable, application of the Baker-Campbell-Hausdorff 
formula gives 

{i’o\A{t)\'ipo) = tT[C^ ad^(A)], (2) 

iV=0 


N times 

where ad:^(A) = [... [A, "H ],...,]?{] is a nested commu¬ 
tator, and Cq is the transpose of the correlation matrix 
in the initial state [Co]ij = (V’o|cJcj|'!/'o)- If the matrix 
Co is known, the representation in power series ([^ can 
be efficiently computed. Due to the locality of A and 
H, only 0{N^) matrix elements in the nested commu¬ 
tator are non-vanishing. One can thus use a finite block 
of the correlation matrix Cq- The coefficients will then 
be exact in the thermodynamic limit as long as the ele¬ 
ments in the nested commutator do not reach the bound¬ 
ary of the block. The method can be applied to de¬ 
rive exact expressions for any two-point function, and 
also for the Loschmidt echo. The correlation matrix in 
the initial state \ijjo) is [Co\ij = C^j for i,j > 0 and 
[Co\ij = for i, j < 0, with matrix elements Cg 

(a = I, r) given by 

=5“(*-j)+5“(*+j), (3) 

and S°'{x) = sin {kpx)/{Trx). The matrix (7“ is a sum of 
a Toeplitz matrix, with entries that depend on i — j, and 
a Hankel matrix, with entries that depend on i + j . 

Long times, semiclassical regime and stationary phase 
approximation - In the limit L —>■ oo the final Hamilto¬ 
nian H = J dk e(k)p{k)f{k) is diagonal in momentum 
space, so the two-point function at time t is 
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FIG. 1. Top. Particle density as a function of x/t for t = 384 
and L = 3072. The filling fractions are kp = tv, kp = Stt/G, 
kp = tt/2 from left to right. Bottom. Comparison between 
semiclassics, Eq. Q, and finite-size numerical simulations for 
a correlator {A„,_At)cx+A{t)). The fillings are kp = Stt/O, kp = 
7r/6. There are three regions. The first \x/t\ < 1/2 is the 
NESS. The correlations in the second region \x/t\ > 1 are 
that of the initial state. Most interesting is the third region 
1/2 < \x/t\ < 1, with position dependent local correlations. 


that it is given by f{k,q) = f\k,q) + f''{k,q), where 
each f/'"(k,q) = frj.''{k,q) + fj^^{k,q). The subscripts 
T and H refer to the Toeplitz-I-Hankel decomposition of 
the matrix These are given by 


fT{k,q) 

fnik, g) 


xijq) + 

27r(l - ’ 

e-^<igi{e-^<i) - C>^gi{A'^) 
27r(e“®^ — e®'?) ’ 


( 5 ) 

( 6 ) 


{clit)cyit)) = I (4) 

where f{k,q) = {tpolPik)f{q)\'<Po) and e{k) = -cosk. 
The integral is taken over [—7r;7r]^. The large-time be¬ 
havior for fixed x and y is then determined m from the 
points where the phase is stationary, as well as possi¬ 
ble singularities in f{k,q). For this type of protocol, it 
is for example known that a NESS develops in the mid¬ 
dle [TUtlb) . Here we are interested in a different regime 
at large time where x/t is kept finite. Such a limit has 
been studied in several particular cases (including the 
DW limit [131 HI] or Fermi gases in the continuum [17l - 
m), but we present a general treatment for all fillings 
here. First we compute exactly f{k, q), the Fourier trans¬ 
form of the initial correlation matrix ([^. We find 


where gi{z) 





. Here x/(fc) = 1 if fc e 


[—kp,kp] and zero otherwise; similar expressions hold 
for f^{k,q). Note that the fj'" have a pole and are not 
analytic at Lkp . With this at hand, the asymptotic be¬ 
havior of Q may be determined by using the stationary 
phase method. The stationary points of the phase in Q 
satisfy the equations v{ks) — x/t = 0 and viols') — y/t = 0 
for general x and y, where v{k) = is the group ve¬ 
locity. In the limit x/t,y/t finite and \x — y\/t \ the 
two stationary points almost coincide and the integrand 
f{k,q) in ([^[^, is singular. The integral is then domi¬ 
nated by the region where k — q is small. Introducing new 
variables K = {k + q)/2 and Q = k — q, the stationary 
point is located at Qs = 0, irrespective of K. Expanding 
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around it and using Jjj ^ q-^o+ ~ ^ix) we obtain 


{cUt)Cy{t)) = 


>-k‘ 


^e-^K(.-v)Q (_x + y 


■ v{K)t 


+ y*' 


( 7 ) 


where 0(x) is the Heaviside step function. The density 
p{x,t) = {c\{t)cx{t)) [ 13 ] is obtained by setting x = y. 
See Fig. top panel; a numerical check at distance 
\x — y\ = 8 is also shown in the bottom panel. As can 
be seen the agreement is excellent and improves when in¬ 
creasing the system size, which confirms that it becomes 
exact in the limit L —> oo. Note that 0 is entirely deter¬ 
mined by the pole contribution in /(fc, g), but subleading 
corrections depend on its full analytic form. 

Such types of results are known as semiclassical (or 
hydrodymamic) approximations in the literature, as each 
fermion at momentum k in the initial state propagates 
ballistically at speed v{k). Eq. Q is then the statement 
that this approximation becomes exact in the scaling 
limit with x/t,y/t finite but \x — y\/t ^ 1. There are 
three regions. In the first {\x/t\ > 1), one of the step 
function is identically zero while the other is identically 
one. The correlations are locally those in the bulk of the 
initial left (or right) ground state. In the central region 
\x/t\ < sin Up the correlations are those of a shifted Fermi 
sea with momenta in [—k’^p^Up], This steady state car¬ 
ries some current. Most interesting are the intermediate 
regions sinfc);. < \x/t\ < 1, which are inhomogeneous. 
In that case the correlations are more complicated, as 
is shown in Fig. [l] Another interesting feature is that 
a naive low energy field theory description would break 
down, due to the inhomogeneous background. Indeed, the 
Fermi velocity depends on position and time, so it has 
to be a field theory in a curved space-time with a met¬ 
ric ds^ = dx"^ — v{x,t)'^dU . Since in l-|-ld all metrics are 
conformally equivalent, it is always possible to come back 
to a flat geometry. However, the relation between lattice 
and continuous geometries is complicated. Such aspects 
are investigated in | 33 |- Semiclassical results can also be 
obtained for other physically relevant initial states, e.g. 
when the two halves are held at different temperatures 
T; and T^. The large x and t behavior is still dominated 
by the pole in /(fc, q) [25] with residue proportional to 
the Fermi-Dirac distributions at temperatures The 
inhomogeneous energy density prohle can be similarly 
computed [22j . A numerical verihcation of the validity of 
a semiclassical picture for energy transport in the con¬ 
tinuum is presented in 

Let us now discuss other large-time regimes. First, it is 
interesting to look at what happens close to the edge of 
the front x/t ~ ±1. For the DW limit it has been shown 
[T^ that suitably rescaled correlations are described by 
the Airy kernel |26j . Such a result follows readily from our 


kp = 7r/3 kp = 2 it/Z 



FIG. 2. Left. Exact expression for the particle current (in 
blue) at fillings kp = it — kp = tt/S obtained from (|^. The 
green curve is the asymptotics expansion (|^ derived with a 
stationary phase approximation. Finally the dashed red line is 
the value attained in the NESS by the current. Right. Similar 
plot for the real part of Q-^(t) = {U_^{t)c 4 {t)) at kp = 2nl3. 


formalism for any filling kp ^ 'k 12 |22|. For large-distance 
correlations with \x — y\/t of order one, 0 breaks down. 
The reason is that the two points k^ and Qs where the 
phase is stationary do not coincide anymore. Our method 
can be generalized to such a case, however, by expanding 
inside the exponentials around this point, and by care- 
ful treatment of the boundary contributions at . In 
this case the computations become more involved, and 
depend on the exact form of f{k,q). Another example 
with boundary contributions is discussed below. 

Corrections to the steady state - We now look at the 
corrections to the steady state predicted by Eq. 0 , in 
the regime xjt I. The simplest example is the current 
Jih) = Im (cj(t)ci(t)) through the link in the middle of 
the junction. From the power series method we obtain 


J(t) = 2^ 


(- 1 )^ 


fc=i 

4fcs — k — s 


(2fc- I)! 


2k-l k r 


E 

S =1 


S\2s) 


( f^k+s — l\2 
\^2k-l ) 


2k- I 


(8) 


with the binomial coefficient (^) and Ju{t) the Bessel 
function of the first kind. Eq. ^ extends to general Up 
the result for the particle current derived in [14] for a DW 
initial state {Up = tt, where only the last term remains). 

The stationary phase treatment of this correlation re¬ 
quires some special care E3, due to the singularities in 
/(fc, q) and the presence of sharp boundaries in A:-space 
when Up ^ TT. The singularity at k = q may be removed 
by studying the time-derivative of 0. The extra fac¬ 
tor e{k) — e{q) gets rid of the denominator in 00 and 
the integrations over k and q are now decoupled [22j . 
The semiclassical results may also be recovered using this 
derivative procedure. The resulting one-dimensional inte¬ 
grals can be evaluated by stationary phase in the large t 
limit and the result integrated back to derive the asymp¬ 
totic expansion of the two-point function. For the current 
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FIG. 3. Left. — ln/!(i) for several fillings a,t V = 0. Right. 
Same with varying interaction strengths V for the DW initial 
state. The length of the system is L = 24. Notice how in the 
gapped phase the Loschmidt echo relaxes to a finite value and 
keeps instead decreasing with time in the gapless phase. 


through the junction we find 


Jit) 


t>i cos{Tr — kp) {tt — 2kp) cos2t 


TT 27r^t 

sin(t_) sin(t cosfcir) + cos{kp — t cos kp) sin(t+) 
2~^Jsm{kp){nt)^J ’ 


(9) 


with t± = t±7r/4. The interpretation goes as follows. The 
approach to the NESS of arbitrary correlation functions 
contains semi-integer powers of t for all kp ^ 0, tt ; these 
fractional powers are due to the aforementioned bound¬ 
ary contributions, and are universal. Also, the frequency 
of oscillations are a; (fcs, A:') = |e(fcs)—e(A:')|, where kg and 
k'g belong to the set of critical points {0, tt , A:^}, obtained 
from the stationary phase. A comparison between the 
exact (|^ and asymptotic @ results is shown in Fig. [2) 

Loschmidt echo - The Loschmidt echo is defined as 
the overlap L{t) = | ('(/'o|e~*^*IV'o) | ■ For short times, a 
power-series expansion gives C{t) « 1 — with 7 = 
(iL^) — (iL)^ the variance of the post-quench Hamiltonian 
in the initial state; 7 can be computed exactly |22j . 

The long-time behavior was studied previously for 
k^F = 7r/2 [2l[3], where only low-energy excitations are 
generated, and conformal field theory techniques (GET) 
may be used. At times much larger than the lattice spac¬ 
ing, it decays as a universal power-law L{t) ^ 
a manifestation of the celebrated Anderson orthogonal¬ 
ity catastrophe [101 m- c is the central charge of the 
GET (here c = 1). In the DW limit, we are able to 
compute C{t) exactly as follows. By application of the 
Wick theorem, C{t) = lim det (gi-j — 5i+,), where 

n—>-oo 

gp = ^ Using a generalization [28] of 

the strong Szegd limit theorem [29] to Toeplitz-I-Hankel 
determinants, we are able to derive 


C{t) = (10) 


This can also be derived using the short-time expansion 
\n\ . This simple formula has several remarkable features. 
First, it is exact at all times in an infinite system (for 


finite L, it is exponentially accurate until t ~ T/2). An¬ 
other is that (the logarithm of) its imaginary time version 
C{iT), T e K, can be interpreted [23] as the free energy 
of a statistical mechanical system where all degrees of 
freedom outside of a circle of radius t/2 are frozen, anal¬ 
ogous to the celebrated “arctic circle theorem” |5D] for 
classical dimers. 

At intermediate fillings, the long-time behavior can 
only be accessed numerically. We find a gaussian decay 
for any ki > 7r/2, i.e., C{t) ~ , with a prefactor well 

described by the formula 7=3: cos^ kp. The short- and 
long-time behaviors of C{t) are summarized in Fig. for 
several fillings, through plots of — ln£(A) versus time. 

Finite interactions - Let us add to the Hamiltonian 
Q the density-density interaction with 

U > 0 and focus on the DW initial state. In the gap¬ 
less phase (0 < U < 1), the stationary state supporting 
a ballistic particle current has been numerically investi¬ 
gated in [m inn 132 ■ Gorrelation functions in the NESS 
are those of the ground-state at half-filling multiplied 
by a space-dependent phase m- In the gapped phase, 
the presence of heavy-mass Hamiltonian eigenstates hav¬ 
ing dominant overlap with the DW initial state [33] pre¬ 
vents the formation of a light-cone and leads eventually 
to absence of particle transport for large times |3I]. All 
these features are also transparent in the Loschmidt echo, 
which we plot in Fig.[^ To avoid recurrences we assume 
t < L/2 and L sufficiently large. In the gapless phase, 
classical free-energy arguments [33] outlined above sug¬ 
gest the large-time behavior — ln£(A) a{V)J, where 
a{V) is a coefficient. In the gapped phase £{t) relaxes at 
large times to a finite L-dependent value. Such a value 
is expected to vanish as L —>■ 00. Relaxation features an 
oscillatory behavior as a consequence again of the pres¬ 
ence of slow bound states with large scalar product with 
our initial state. 

For small times, the Loschmidt echo behaves as 1—A^/4 
for all V [331 |3Sj , since the only process contributing to 
the cumulant {H'^) — (i^)^ corresponds to moving the 
rightmost particle one step to the right, and then back. 

Conclusion- In this Letter we presented a unified for¬ 
malism to characterize both the short- and long- time 
limit of correlation functions in a fermionic chain. Our 
main example was a system with a inhomogeneous den¬ 
sity profile but the method applies to a wide class of 
initial states. Together with exact expressions for the 
particle current and the Loschmidt echo, we provided 
a full analytical derivation of the so-called semiclassical 
regime and outlined its domain of validity. We showed 
how subleading corrections to correlations at the edges 
of the front are described by the Airy kernel, giving sup¬ 
port to a connection between large deviation functions for 
transport problems in fermionic systems and the Tracy- 
Widom distribution. Our work raises a number of inter¬ 
esting questions: among them, the possibility of a field- 
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theory formulation of the scaling behavior, the study of 
the entanglement spreading in a inhomogeneous back¬ 
ground and a deeper understanding of interaction effects, 
which might be achieved by engineering a tailored Bethe 
rapidity distribution for the NESS. 
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This supplementary material contains additional information about 

• The calculation of the function f{k,q) in the inhomogeneous density quench (Sec. 0 

• The power series expansion of correlation functions (Sec. H 

• The stationary phase approximation of correlation functions (Sec. 0 

• The Loschmidt echo (Sec. 


A. CALCULATION OF THE FUNCTION f{k,q) IN THE INHOMOGENEOUS DENSITY QUENCH 

Before connecting them the two fermionic Hamiltonians have open boundary conditions. In the thermodynamic 
limit L —^ oo, they are separately diagonalized in momentum space by the following transformations 


Cj = J dp sin{pj)fr{p) j > 0, 

Cj = J dp sm{p{j - l))fi{p) j < 0, 


(51) 

(52) 


where the fermionic operators /^^j(p) create a fermion with momentum p on the right and left half of the chain 

respectively. The one particles states before the quench are therefore \p)r/i = fl/i{p)\0)- After the quench, the Hamil¬ 
tonian is fully translation invariant and diagonalized by Fourier transform, i.e. the local fermions Cj are expressed by 


^ /A dk e^^^ f{k) and the one-particle state after the quench are given by \k) = p{k)\0). The matrix elements 
between the single-particle states before and after the quench are 




Mr{k,p)= 

Mi{k,p) = {k\p)i = ^ 

ZTTl 


1 




-{p-^ -p) 


o-W 




- (p-)- -p) 


(53) 

(54) 


The initial state |'0o) we chose is the factorized Fermi sea \'ifi)\'tpf) and one simply has (V'ol/l(p)//3(E)ll/’o) = S{p— 
p')5a,pQ{kf, - p), for a, 13 = {r, 1}. Then f{k, q) in Q is the sum f{k, q) + f''{k, q) where 


f°‘{k,q)= dp M*{k,p)Ma{q,p). 


(S5) 


The integration domain in (S51 can be extended to p G [~kp, kp ], exploiting the symmetry properties of the functions 
in (S3 S4). It is also useful to further decompose fa{k,q) into the sum fp{k,q) +/))(A;,< 7 ). For example one has 


fxik, q) — 

fnik^q) = - 



I 

(S6) 

(1 — _ gi(9-p+*o)) ’ 


g2ip 

(S7) 

(1 — g-i(fc-p-*0)^q _ gj(9+P+*0)) ’ 


and very similar expressions for f'"{k, q). Notice that the subscripts T or H refer to having chosen products of terms in 
(S3 S4), singular for k = ±p and q = ±p (subscript T) or k = ±p and q = =FP (subscript H). The integrals in ( S6|S7 ) 
can be computed changing variable z = and performing integration along the contour A depicted in Fig. SI Only 
the singularities inside the contour contribute to the hnal result. 
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FIG. SI. Integration contour A utilized to evaluate the integrals in (S6 S71. Consider for example the function frik, q) in (S 81 , 
then the only singularity inside the contour A is aX p = q -\- iQ. 


Performing those integrations we obtain 


dTT fT{k,q) = 


_ gi{q—k+iq) 


log 






- log 


gikp _ ^iq 
_ gjg 


- 27rfxi(g) 


g-.^ _ g^g 

I 

^ _ Qi{q—k—iO) 


'log 


-e^Mog 


g '^kp g —iq 




log — 




g "^kpi _ 


pOn.p - „ 

— log I - T— -- 

^ \ g—2/c^ _ ^tq 


2'KiXrik) 


4TT'^f^{k,q) = — 


ik 


- 


^ — ik 


log — 


^ik'^F _ ^—'i'k 


g ^kF _ g —zk 


- -e^nog — 


g ^kF _ 


(58) 

(59) 

(510) 

(511) 


where Xa{x), is a function with value one if a; € [—and zero otherwise. The initial state correlation matrix 
[C^o] nm ii’ol IV'o) is given by 


[Co\nm. = Q{n — 1)0(to — l)[5''^(n — m) + S'"{n + m)] + S{—n)Q{—m)[S\n — m) + S\—n — m + 2)], (S12) 

where S°{x) has been defined in the main text below It is possible to further check ( S^Sll ), showing that they 
yield indeed the Fourier transform of [Co]nm, 


[Co]nm = ^fdk y" dq 

as it follows from Q. To verify (S13) it might be useful the following identity valid for s S Z and a G 

sin(as) 


dz z""Mog — 


e — z 


— z 


= 47r0(s — 1)- 


(S13) 


(S14) 


with C the unit circle. 


B. POWER SERIES EXPANSION OF CORRELATION FUNCTIONS 

An exact expression for the coefficients in the power series expansion of the correlation functions in ([^ can be 
obtained starting from the knowledge of the nested commutator ad^(A). For the Hamiltonian Q and A = c|.Cy 
simple combinatorics leads to 


[aA^{A)]ij 


(- 1 ) 




2N 


N 

N+\d-{x-y)\ 

" 2 " 


N 

N+\D-{x+y)\ 

. ^ „ 


(S15) 







































kf-p = '^/3 



FIG. S2. Real part of the two-point functions Qx{t) = (cl 2 ,(t)ca;-i-i(t)) plotted for different values of |x| = 2, 5 at kp = n — kp = 
tv/3. The curves are obtained from the power series expansion ( |S17[ ) summing N = 255 coefficients. We observe clear signatures 
of light-cone effects at short times and relaxation to the stationary value in (S19l, denoted by dashed lines in the figure, at 
larger times. 


where d = i—j and D = i+j, satisfying the constraints \d—{x — y)\ < N, \D—{x+y)\ < N and N—{x — y) = D mod 2. 
Moreover when y = —x -I- 1 (a; < 0) the points are taken symmetrically on the left and right halves and it is not 
difficult to prove that all the coefficients in § are vanishing if < |a:|, which is a consequence of light-cone effects. 
From ([^ it is clear that 

{cl{t)cy{t)) = + an, (S16) 

N=0 


where 
case X = - 


= Tr[C'J ad^(M)] and the matrices Ca 
-y -\-l and |a;| = —x\ then one has 


defined in (S12) for a 


l,r. Let us focus again on the symmetric 


Af-2|a:|-l 




= = T! (- W ( 7V+M^|.| + 1| ) (" 1) ( iV+D-1 ) 

n,m>l d=-N+l ^ 2 ' D=2+\d\ V 2 / 

AT-HI min(n-2.7V-2|a;|-l) / AT \ 

+ (-!)'"' E' E' 

n—o_Lvv^-^^/o\ N 2 / _/'n_o^ V 2 / 


(S17) 


D=2+mod(Af-|-1.2) 


d=-{D-2) 


with the primed sum restricted to d = iV -f 1 mod 2 and D = N + 1 mod 2. The expression (S17) can be simplified 
in the case a; = 0. The coefficients are the coefficients in the power series expansion of the particle current J^ 

we get 


fc-i 


a. 


(2fc-l) _ 


= -E 


s=0 


(fc-bs)5''(2s) 
2fc- 1 




fc+s \2 


(S18) 


where the primed summation denotes now that the term containing Sr(0) has to be divided by two. We also introduced 
the notation for the binomial coefficient (^). Analogously, one hnds the odd coefficients and summing the 

two, see (S16), we derive the result given in ® for the particle current. Although (S17) appears complicated the 


coefficients of the power series expansion can be derived in symbolic form for all two-point functions. We show in Fig. 
S2 typical plots of the symmetric correlators Gx = (^o|cl;(t)c_a;+i(t)|^/’o) obtained with the method above. The signal 
starts to be non-vanishing for t ^ |a;| which is a clear indication of light-cone effects. Correlation functions approach 


their stationary value at long times (dashed in curves in Fig. S2) 


(c^c - ^-^su-y) Mko{x - y)) 

{CxCy)NESS 


(S19) 


for 6 = kp — kg and ko = tt 12 that will be derived in Sec. C2. 
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C. STATIONARY PHASE APPROXIMATION OF CORRELATION FUNCTIONS 


The asymptotics of correlation functions {cl;{t)cy{t)) can be studied mainly in three interesting regimes determined 
by the value of the parameters x, y and t entering the Fourier transform Q. We will consider the cases 


x/t ^ 0 and y/t ^ 0 with t ^ 1, see Sec. Cl The points x and y are inside the light-cone and correlations 


have already relaxed to their stationary value in the translation invariant NESS. The approach to the stationary 
state is oscillatory with corrections organized in power series of 


x/t and y/t finite but all the variables large, see Sec. C2 This particular limit is known in the literature as the 


semiclassical (or scaling) limit; we offer an alternative derivation of our main results in Sec. C3 


x/t ~ ±1, see Sec. C4 We focus on the behavior of correlations at the edge of the fronts and show that they 
are described by the Airy kernel. 


Cl. Corrections to the steady state 


We start considering the first case that technically amounts to performing a stationary phase approximation in the 
variable t of the double integral Q. The function /(fc, q) is singular for k = q and we found it easier to consider first 
the time derivative of the two-poin t correla tion function. The extra factor e{k) — e{q) exactly cancels the divergence 
in both f!(/'~{k, q) and fJi^{k, q) in ( S^Sll ). Moreover the integrations over k and q are now decoupled and one needs 
to analyze the large-time behavior of the following one-dimensional integrals 


A{n,t) = — dke 

27r 7-,r 

f'kp 

f-kp 


it cos k ink 


1 f ^ 

B{kp,n,t) = / dke 

J-kp 

1 r 

C{kF,n,t) = TT- / dke 
J-p 


it cos k ink 


it cos k Ank 


e“«log — 


^ikp _ 


-ikp 


oik 


(520) 

(521) 

(522) 


which is a more comfortable situation. We will provide asymptotic expansion of correlators up to 0{t~^^^), that implies 
that we have to keep terms in the asymptotics expansions of the derivative up to 0(t~^^^). Since the co ntributio n of 
ordinary stationary points is of the form for p positive integer, we have to expand the functions ( S20|S22 ) up 

to order The first integral (S201 is proportional to a Bessel function whose asymptotics are known. However it is 


good exercise to recover it. There are two stationary points at fc = 0 and k = n and the result is 

e-/2y2sin(f-^+t) 4n^_ie-/2y2sin(f+ ^-t) 


(S23) 


Notice that the asymptotics expansion of A{n, t) only contains seminteger powers of t. To derive exact asymptotics 
for the function B(kF,n,t) one has to add the boundary terms at fc = ikp to the contribution from the stationary 
point located at fc = 0. Boundary contributions can be obtained by making the change of variables u = cos k and then 
integrating by parts m- These boundary terms generate integer powers of t in the expansion 

cos nkp — 1 (nsinnkp + cot kp cosnkp) r/,. 

B{kp, n, t) = —-b-——-A +-^--^- 0 - 2 , - - - — + 0{t-^'^). 


\/Ari 


Trt sin kp 


8 y/ArfiA 


sin kp 


(S24) 

Then we are left with the problem of finding the asymptotics expansion of C{kp, n, t). It is convenient to split it into 
two parts C = Cl + C 2 where Ci and C 2 contain the imaginary and real part of the logarithm respectively. We have 


Ci{kp, n, t) = ikpA(n, t) — ipB(kp^ n, t) 


and the asymptotic follows from (S23)-(S24), whereas 

1 


C2{kF,n,t) = 


47r 


dk git cos k^^nk 


log 


1 — COs[fcF — k] 
1 — COs[fcF + k] 


(S25) 


(S26) 
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The asymptotics is dominated by two ordinary stationary points at /cg = 0, tt where the phase vanishes; they produce 
terms since the integrand is vanishing there. To understand the contribution of the logarithmic singularity 

at k = ±kF we rewrite the integral as 


C2{kF,n,t) = — 


dfce**“"''sin(nA:) log 


1 — COs[fc F — k] 


(S27) 


;0 V 1 ~ COs[fci? + k] _ 

and notice that we can expand the logarithm for fc —as 21og(fc — /cf) + reg.. When crossing the logarithmic 
singularity the imaginary part jumps by 2z7r and generates an effective boundary contribution at that can be 
evaluated by integration by parts. One finds the final result 


C 2 {kF, n, t) =e 


it cos kp f i sin nkF coikF sinnkF — ncosnkF 


t sin kF 


+ 


sin"^ kF 


An+l 


TT (sin L 


/TT riTT 

coskf (^4 ^2-^ 


^ f sin ^ 


TT riTT 

- \-t 

4 2 




(S28) 


The asymptotics expansions ( |S23[ ), ( |S24[ ) and (S28) solve our problem as the derivative of any correlator can be 
expressed in terms of them. We decompose the two-point function into the sum of four pieces r, I (right/left) and T, H 


(Toeplitz/Hankel) according to the function f!Jljj{k,q) entering in each of them, see (S8 


Sll), and get 


dt{ci{t)cm{t))H,i = - — [-{A{-n,-t)C{kp, 1 - m,t) - A{fn,t)C{kp, 1 - n, -tj) 

-1 47r 


-|-A(1 — n, —t)C{kF, 2 — m,t) — A{m — 1, t)C{kF, 2 — n, —t)], 


(S29) 


dt{cl{t)cjn{t))T,i = -^[Cikp, -n + 1, -t)A{m,t) - A{-n + 1, -t)C{kp,m,t) - C{kp, -n, -t)A{m - l,t) 

_^_-_I 47 r 


-I- A{—n, —t)C{kp, m — l,t) — 2'Ki{A{—n -\- 1, —t)B{kp, m, t) — A{—n^ —t)B{kp,m — 1, <))], 


(S30) 


dt{cl(t)cm(t))H,r = “ {A{m, t) C{kp, l + n, -t) - A{-n,-t)C{kp, 1 -b m,t))+ 

-1 47r 


A{m — 1, t)C{kp,n, —t) — A{—n -b 1, —t)C{kp, m, t)], 


(S31) 


dt{c\^{t)cm{t))T,r =^[C{k''p, -Ti -b 1, -t)A{m,t) - A{-n + 1, -t)C{kp,m,t) - C{kp, -n,-t)A{m - l,t) 
--—--1 47r 


-b A{—n, —t)C{kp, m — l,t) + 2Fi{B{kp, —n + 1, —t)A{m, t) — B{kp, —n, —t)A{m — l,t))]. 


(S32) 


Plugging (S23), (S24| and (S28) into ( S29|S30 ) we determine the full asymptotics expansion of the time-derivative of 
the two-point functions. The result can be integrated back to derive the approach to the NESS of all the correlators. 
We gave a relevant example for n = 0 and m = 1, in the main text Eq. ([^. To validate the method we also provide 
the asymptotic expansion of the real part of the correlator G-sit) at = tt — = 27r/3. It reads 


un ^ , -8/v^ -b TT 7 cos(2t) 
n5-3it)] c. _ + - 


2cos[g(7r — 9t)] — 3(3 -b -s/S) cos(t/2) -b cos(3t/2) -b [—9 -b 2-\/3 — 2-\/3cos(t)] sin(t/2) 

6(7rt)3/^ 


. (S33) 


The comparison between the asymptotic formula (S33l and the exact curve obtained by the power series expansion, 


discussed in the previous section, is shown in the main text. 


C2. Semiclassical limit 


We now obtain results for the correlators in the semiclassical limit. The stationary points of the phases in Q of 
the main text are the solutions of the following equations (u(fc) = de{k)/dk) 

v{ks) — x/t = 0, (S34) 

- y/t = 0, (S35) 
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X 


FIG. S3. Semiclassical limit for the energy density profile h{x,t) plotted at different times t and for f3i = 1, Pr = 0.5. The 
function ( |S41[ | is expected to exactly reproduce the large x, t limit of the energy density in a fermionic chain, whose halves are 
at time t = 0 thermalized independently at temperatures = Pi/l-- 


We restrict ourselves to the case {x — y)jt ^ 0; taking the difference of ( S34|S35 ) we observe that stationary points 
either coincide = kg or are shifted hy tt, kg = t: — qg. Let us introduce new variables Q, K, d, X through 


k = K + Ql2, q = K-Ql2, 
x = X-dl2, y = X + d/2. 


(536) 

(537) 


When the saddle points coalesce the integrand in Q is singular and linearization around Qg = kg — qg = 0 gives the 
leading behavior in the semiclassical region. We obtain 


{cUt)Cy{t)) 


semiclassics 


27r 


dK 


-k‘ 


ei-x + viK)t) + — 


dK 


Q{X-v{K)t), 


(S38) 


-kl 


where we have used the standard representation of the Heaviside step function 0(x) = ^ Ir^Q q^^q - For X/t ^ 0 
we are back with correlation functions in the translation invariant NESS. After elementary integration one finds 


(4cy)NESS = 


I ^ikp(x-y) _ ^-ik‘p(x-y) 

2Tri X — y 


(S39) 


the particle momenta in the NESS are obtained shifting all particle momenta in the ground state k G [—fco, fcg] by a 
constant amount 6 = kp — ko (ko = 7r/2), see (S19). 

As we discussed briefly in the main text, semiclassical results can be obtained in other physically relevant contexts. 
For example, when the two halves of the chain have same densities but different inverse temperatures /?; and /3r one 
expects, after the quench, a ballistic energy transport with a inhomogeneous energy density profile h{x,t). In the 
semiclassical limit, the function h{x,t) is again a scaling function h{x/t) that can be determined analytically. Let po 
be the initial density matrix for the two disconnected chains then [25] 


f{k,q) = TT[poc\k)c{q)] = - + regular 


47r(l _ gi{q-k-iO)'j — gi(q-k+M)'^ 


(S40) 


where we have omitted regular terms for fc —>■ g and ni/^{k) = 1/[1 + Notice that for cosine dispersion 

relation the choice of a zero chemical potential fixes the particle density to be homogeneous and equal 1/2, therefore 
particle transport is absent. The similarity between |S40|) and (pSl, dSTol) is evident. Proceeding in complete analogy 


with the calculation of the correlation functions, explained few lines before, we obtain for the energy density profile 


h{x, t) 


semiclassics 


1 


1 


— / dK e{K)ni{K)Q{—x + v{K)t) + — / dK e{K)nr{K)Q(x — v{K)t) 


2tt 


27r 


(S41) 


where we assumed e{k) = — cos k. A plot is given in Fig. S3 Nothing prevents a study of the large-deviation function of 
the energy-flow, for example with /?; = 0 and /?r = oo, in the same spirit of |12j . Those aspects will be not investigated 
here. 
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C3. An alternative derivation of the semiclassical limit 

We provide here another derivation of our main result, that is closer to the standard stationary phase procedure. 
As already emphasized the main complication is that the stationary points in k and q almost coincide, and the 
denominator is singular in such a limit. This difficulty may be circumvented by using the same derivative trick as in 
section [cn 

Let us explain the method on the term involving flp{k,q). Since f\j{k,q) is regular, it will not contribute in the 
scaling limit. For the same reason, only the pole in /y(fc, g) matters. Therefore, to the leading order we have to 
evaluate 


4WsW)z = 


dk 


dq e 


[£{k) — £{q)]t—ikx-\-iqy 


l-T, 2 tt 2tt 1 - eh9-fe+io+) 

Now we consider the time-derivative of the previous equation. We obtain 


(S42) 




x\ y\ J/i 2 ty 27r 1 - ed9-fc) 


(S43) 


which is now regular at k = q. The phase ^xik) is given by 

^xik) = e{k)t — kx. (S44) 

We now have to look for the points where the phase is stationary. These are solution of the equation 

v{ks) = j (S45) 


where v{k) = = sinfc. For \x/t\ < 1 there are two real solutions 


kt = arcsin — 
t 

The treatment for $j,(g) is similar, and we get 


+ ■ y 

qj = arcsin - 
t 


kg = TT sign X — arcsin —. 


_ . . y 

qg = TT sign y — arcsin -. 


(S46) 


(S47) 


Let us now assume that the kf , qf belong to the integration domain, which means we are in the inhomogeneous 
region —1 < xjt^yjt < — sinfc^ Around the stationary points = 0, the phase may be approximated by 


= ^x{kt) + - kfr, 


(S48) 


and the same goes for ^y{q). There are in principle four contributions at (fc+, g+), {kj,qj), {kf, qj), (fcj, g+). One can 
check that the last two give fast oscillating contributions, which become subleading when integrated back. Neglecting 
those two we obtain 




+ {kt ^ kj ; g+ ^ q^}. (S49) 


/4>"(fc+) 

At distance \x — y\/t <C 1 we have fc+ = qf (resp. kj = qj) to the leading order. Therefore, 






sin(fc+) + e^(v-^)K 


2'Kt cos(fcj) 


(S50) 


(S51) 
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Noticing that tan/c^ = , this may be rewritten as 


d , , , ,, dkt dk- e^<^y-^)K 


dt 


rkt 


dK 


-iK(x-y) 


i.- 2 tt 


(552) 

(553) 


which is the expected result in the left inhomogeneous (front) region —1 < xjt,yjt < —sinkp, see Fig. S4 Note 


however that the method only gives the correlations up to some integration constant, which depends on x and y. 
This constant may be fixed by using a continuity argument. Indeed by applying the same method, we find that the 
derivative vanishes in the regions \x/t\ > 1 and \x/t\ < sinfc^. Since we know that far on the left (resp. right) the 
initial correlations are that of t he in itial left (resp. right) ground state, we use the boundary condition at t = |a;| to 
fix the integration constant in (S53). We then do the same at \x/t\ = sinkp to obtain the correlations in the region 
\x/t\ < sinfc^. 



FIG. S4. The picture shows the integration domain in the variable K constrained by the two Heaviside step functions 0{—x + 
v{K)t) and 0(x — v{K)t) in ( |S38[ ), for x = y and fixed time t. The black curve K{x) is the solution of the stationary phase 
equation v{K) = xjt. The stationary points (Eq. ( |S46| ) in the inhomogeneous front region are shown in blue. Approaching 
the right edge of the front kp — 'k/2, the turning point. Near the turning point, i.e. for x/t ~ 1, correlations 

have non-trivial subleading corrections with respect to their value in the bulk of the right ground state; these corrections are 
described by the Airy kernel (S56l. 


C4. Correlations near the edge and the Airy kernel 


Let us focus for simplicity on the the right edge of the front and assume kp > 7r/2. At the leading order correlations 
are given by their ground-state value in the right-part of the chain and are in particular time-independent. However, 
approaching the boundary of the light-cone the stationary points and qf in ( S46|S47 ) both reach the turning 
point kp = 7r/2 (see Fig. S4) where the second derivative of the phases in (S49l vanishes. As a consequence one finds 


non-trivial subleading corrections that can be computed evaluating the integral 


{cl{t)Cy{t))l = / — 


dk 


-fc 


27r 


\ _ gi{q—k+iO) 


when k,q ^ kp. Expanding up to third-order in the phase, we get 

roo jZ coo j ;: .^ — ik(x — t) — ii^+iq{y—t)+ii^ 


cl{t)Cy{t)) 


near the^ edge — 


dk 


dq e 


I —oo J —oo 


i{k — q — iO) 


(S54) 


(S55) 
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with k = k — and q = q — k^. Introducing the scaling variables X = {x — t) (f) Y = {y — t) (|) and defining 
K = {t/2y/^k, Q = (t/2)^/^q we can rewrite (S55) as 


where K{X,Y) is the celebrated Airy kernel [55] 


-iKX-iK^/3+iQY+iQ^ 13 


27r i{K-Q-iO) 


KiX,Y) = 


Ai(A)Ai'(r) - Ai'(A)Ai(r) 


X -Y 


(S56) 


(S57) 


The last passage in ( |S56| ) follows from the relation —(9x +9y)A'(A, Y) = Ai(A)Ai(y) and the integral representation 
of the Airy function Ai(A) = ^^tXK+iK /3_ Extension of the result to finite temperature correlation functions 
(see the end of Sec. C2) is also possible. 
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D. LOSCHMIDT ECHO 


This section is devoted to the study of the Loschmidt echo after the quench from the domain-wall initial state. We 
provide two separate derivations (Sec. D2|D3 ) of the result 

C{t) = (S58) 


quoted in the main text, as well as a generalization to any dispersion relation. Before doing so, we first briefly discuss 
the short- and long-time behavior of the Loschmidt echo for arbitrary hllings. 


Dl. Short- and long- time behavior 


For convenience we start from the definition C{t) = |('!/'o|e Pi where Hq is the pre-quench Hamiltonian. 


If |'0o) is an eigenstate of Hq, this is clearly immaterial. By expanding for small t one obtains C{t) = 1 — 
with 


7 = Var(iL - Hq) = - i?o)'IV'o) - [ii’oliH - iLo)|V'o)]'. (S59) 

If we have a quadratic operator A which is represented by the matrix A in the fermion real Fock space then the 
variance of such operator in the initial state j'i/'o) is 

Var(H) = Tx[C^ACA], (S60) 

where C = \ — C and C is the initial state correlation matrix. For the filling fraction quench the matrix C is given in 
the main text Eq. ([^. It is possible to check that the values of 7 are independent from the cut-off L, used to represent 
the matrices H, Hq and C and we find the exact expression 

. 27r^ -I--b 27rsin2a: -b sin^(2a;) — 4a;(7r -I- sin2x) 

IA) = - — -, (S6I) 


where x = kp. A numerical check of formula (SGI) is presented in Fig. S5 


As already commented in the main text the Loschmidt echo large time behavior is unfortunately accessible only 


numerically. The data are compatible with a gaussia n de cay C{t) 
empirical formula jikp) = \ cos^(fc^), see again Fig. S5 


_ 


e , with coefficient 7 well reproduced by the 




kp/u 


kp/u 


FIG. S5. Left. We show numerical data for the coefficient 7 (red dots) compared with the exact formula (S6I1. Right. The 
red dots show the coefficient 7 for the farge-time fimit of the Loschmidt echo extracted by the numerics. The bfack fine is the 
empirical formula ^{x) = 4 cos^(a:7r). 


D2. Domain wall limit: determinant derivation 

The first strategy is to consider a finite system of size L, and to express the Loschmidt echo as a determinant. This 
determinant is then evaluated in the limit L —>■ 00. The Hamiltonian can be put in diagonal form 

= (S62) 

q 
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with 


dl 




±s.n 

1 = 1 


t 

L+ 


(S63) 


Here g G {1, 2,..., L} 
to shift the labels of 
convention the initial 

by 


and Eg = cos More general dispersion relations will be discussed later. Note that we chose 
the sites by L/2, so that the open chain now starts at site 1 and ends at site L. With our 
state is IV'o) = Il^ii |0). In a Heisenberg picture, the real space fermion operators are given 





L + 1 


’*4 


2 

L + 1 


L / L 

E E- 


L + 1 


sin 


L+1 



(S64) 


In the limit L —>■ oo for fixed t, the sum over q in the previous equation converges to a sum of two Bessel functions, 
so that 


CfW = E (b65) 

where Jx(t) = cosq-iqx jg ^ Bessel function of the first kind. Note the appearance of Jj+i{t), which comes 

from the boundary. Using Wick’s theorem, the echo may be expressed as a determinant. In the limit L —i oo, we 
obtain 




lim 

L—^oo 


det - Jj+i{t)) 

i<3,l<L/2 


(S66) 


It is important to understand that we have first taken the limit L —>■ oo to get the Bessel functions from (S64), and 


only then taken the infinite determinant limit. Hence the calculation is not mathematically rigorous. The physical 
justification of this assumption is the presence of a light-cone, which ensures for any finite time t that the effects of 
the boundaries at distance L/2 ^ oo are suppressed. 

The r.h.s of the previous equation is the determinant of a matrix whose elements only depends on j — ^ and j + 1. 
Such determinants are called Toeplitz (j — 1) + Hankel (j + 1) determinants. They have been widely studied in the 
mathematical literature, starting from the work of Szegd [Slj . The important object to consider for such evaluations 
is the symbol of the determinant, namely the inverse Fourier transform of Jx(t). Here the symbol simply follows from 
the definition of the Bessel function, and is 


gik) = Jx{t)e 
x^lL 


—ikx _ cos k 


k G [—TT, tt]. 


(S67) 


We note [g{k)]x = = Jx{t) the corresponding Fourier coefficients. With this at hand, the asymptotics 

of (S66) follow from a theorem derived in Ref. [S2]. For sufficiently regular symbols, the following asymptotic formula 
holds 


lim det i[gik)]j-i - [gik)]j+i) = exp - V x {[\ogg{k)]xf - V[logg]2a; 

n—>-oo l<j,l<n \ Z 


(S68) 


provided g{—k) = g{k) and [log5(A:)]o = 0, which is the case here. The (logarithm of the) symbol (S67) corresponding 
to the Bessel function has only one non zero Fourier coefficient, so that we get 






(S69) 


which is exact at all times in an infinite system. Squaring this gives the claimed result (S58) for the Loschmidt echo. 


It is important to emphasize that the choice of boundary conditions is not innocent in such calculations. In case of 
periodic boundary conditions, there are two light cones developing after the quench, one at x = T/2, L/2 + 1, but also 
one at X = 0, L. Because of these two light cones, the Lochmidt echo in a periodic system is the square of the one in 
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an open system. This observation serves as an alternative and simpler way to compute it. Indeed with PBC we get a 
Toeplitz determinant 


n—^oo l<j,l<n 


= exp E a;[logg(fc)]^[logg(fc)]_^ 




= 

AtH. 




(570) 

(571) 

(572) 

(573) 


To go from (S70) to (S71) we have use the strong Szego limit theorem [ST]. Generalization to other dispersion relations 


is also possible. Let us consider an Hamiltonian with longer-range hoppings 

p 




(S74) 


The range p of the longest hopping term need not be finite, but let us assume that for now. Then, regularizing the 
Hamiltonian on a periodic ring with L sites, the dispersion relation is 


5(fc) = 2^ Ua cos{ak)^ 




s(afe) 


so that the symbol of the Toeplitz determinant becomes 

Application of the Szego limit theorem then gives 


Taking the square-root we obtain 


and so 


(e*‘^)=exp 


C{t) = exp ( — 


E' 


E' 


(575) 

(576) 

(577) 

(578) 

(579) 


We have assumed that p is finite in the calculation, but the result should also hold for an infinite number of Fourier 
coefficients, provided the series X]a>i converges. 


D3. Domain wall limit: combinatorial derivation 

We provide an alternative derivation that uses only combinatorial means. Let us focus on the dispersion relation 
e(fc) = ucosk for now. The method is similar in spirit to the short-time expansion presented in the main text, and 
can be made fully rigorous. The real-space Hamiltonian is 

“ E + c]g+i) ■ (S 80 ) 

jez 

Expanding in power series yields 

= (S 81 ) 

m—0 
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where (.) denotes the average in the domain wall initial state. Therefore, evaluating the Loschmidt echo boils down to 
finding an exact expression for all the moments (iJ™) of the Hamiltonian. The initial state can be pictured as shown 
below: 


- - •-•-•-•-•- m- 


-5 -4 -3 -2 -1 0 


1 


2 


3 


4 


5 


« - 
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The sites on the left (resp. right) are filled (resp. empty). The action of H on \tpo) = na;<o 1^) obeys rather simple 
rules, which we explain now. Each particle tries to go to one of it’s nearest neighbor site, provided there is not already 
a particle here. At the first time step the only possibility is for the rightmost particle to go one step to the right. 
Hence 


H I'tPo) =H\... 111111|000000...) = w |... 111110|100000...). 


(S82) 


Here a 1 represents a site occupied by a fermion, and a 0 an empty site (hole). To improve readability a vertical bar 
I between sites 0 and 1 was added. is given by 

l^g) =Hu\... 111110|100000...) (S83) 

= u^\... 111111|000000 ...)+u^\... 111110|010000 ...)+u^\... 111101|10000...) . (S84) 

Since {H'^) = ('0o|^^”^|'0o) selects the initial state in this expansion, we have (H) = 0 and (i?^) = u^. Said differently, 
(H™) counts the number of ways for hardcore particles to move to the left and to the right, where at each step only 
one particle moves, with the constraint that all particles have to come back to their initial positions after m steps. 
Because a particle is forced to move at each step, it is easy to see that = 0 for integer m. With these rules 

at hand, obtaining boils down to a known combinatorial problem. The result is 

(S85) 


where (2m — 1)!! = 1.3.5... (2m — 1) is the double factorial. This formula is remarkably simple, but its derivation is 
not. Proving it can be done by (i) reformulating the problem in terms of oscillating Young tableaux, and then (ii) 
using a known bijection between those and perfect matchings. 

The first step (i) goes as follows [S3]. To each particle configuration may be associated a Young tableau, i.e. a 
collection of left-justified boxes, with non increasing row height. This is done by first removing all ones part of a 
semi-infinite sequence of ones, and doing the same for the zeroes. The height of the first row is then the number of 
remaining ones (which equals the number of remaining zeroes). We then draw an horizontal edge to the right (resp. 
vertical edge up) for each zero (resp. one) encountered. The initial state I'l/'o) corresponds to an empty diagram. This 
procedure is illustrated in Fig. S6 on the example {ip) = |... 1111100101|1100100000...). Moving a particle then 


corresponds to adding or removing a unit box. 





1 




^0 0 
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1 ° 


0 0 



FIG. S6. Young tableau corresponding to the state \ip) = |... 1111100101|1100100000 ...). The relevant particle configuration 
after removal of the spurious ones and zeroes is 00101| 11001. The initial state corresponds to the empty diagram. 

Evaluating (H^™) amounts to counting sequences of young tableaux that start and end at the empty diagram, with 
m steps up and m steps down. Such sequences are called oscillating tableaux, and they have been shown to be in 
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bijection with perfect matchings |S4j . Since the number of perfect matchings of (1, 2,, 2m) is exactly {2m — 1)!!, 
the result of Ref. [S4] immediately implies (S85). Hence 


{«'"> = 1 + E 

m—1 


(2m)! 


(2m 


1 )!! 


Squaring this and setting u = 1/2, we finally obtain 


L{t) = 




(586) 

(587) 


(S88) 


as claimed in the main text. The technique can also be generalized to other dispersion relations. In that case the longer- 
range movements of the particles correspond to the addition of bigger boxes, called ribbons, with signed weights that 
depend on the number of particles overtaken at each step. This exactly encodes the free fermionic nature of the 
problem. Such oscillating ribbon tableaux have been studied in Ref. [S5| by similar bijection techniques. Using the 
bijection described in the reference, it then follows that |S6] 

^ (2m - 1)!! [ul + 2ul + + ...] , (S89) 

which implies 


C{t) = exp 



(S90) 


This result coincides with the determinant method. Note that the Loschmidt echo is zero at all times t > 0 in case 
the sum X]ct>i diverges. 
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